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Obtaining general relations between macroscopic properties of random assemblies, such as den¬ 
sity, and the microscopic properties of their constituent particles, such as shape, is a foundational 
challenge in the study of amorphous materials. By leveraging existing understanding of the random 
packing of spherical particles, we estimate the random packing density for all sufficiently spherical 
shapes. Our method uses the ensemble of random packing configurations of spheres as a reference 
point for a perturbative calculation, which we carry to linear order in the deformation. A fully 
analytic calculation shows that all sufficiently spherical shapes pack more densely than spheres. 
Additionally, we use simulation data for spheres to calculate numerical estimates for nonspherical 
particles and compare these estimates to simulations. 


Understanding the relationship between macroscopic 
properties of a random assembly and the microscopic 
properties, such as shape and composition, of the par¬ 
ticles that make it up is a fundamental problem in ma¬ 
terial science. Colloidal suspensions [ 1 ], granular ma¬ 
terials [2], nanoparticle assemblies [3], and biomaterials 
mg all feature disordered agglomerations of nonover¬ 
lapping particles. When the particles are all undeformed 
spheres of the same size, the typical maximal density, 
measured by volume fraction, obtained by random pack¬ 
ing is (j) ~ 0.64. Though the precise value of the density 
depends on preparation protocol, the conventional value, 
to within a percent, has been reproduced in many exper¬ 
imental and computational settings Hi]- This robust 
value suggests a purely geometric phenomenon, irrespec¬ 
tive of particular system-specific details such as residual 
interaction, hydrodynamics, or gravity, and this density 
is known as the random-close-packing density. The term 
random packing implies the lack of crystalline order in 
the packing, even at short ranges. Without this restric¬ 
tion, the largest volume fraction obtained by packing of 
nonoverlapping spheres is 7 r/ 3\/2 « 0.74 [TO]. 

While the case of equal-sized spheres is the canoni¬ 
cal one, applications in granular materials mm, bio¬ 
physics m m, liquid crystals m, and self-assembly 
[niini have driven interest in the random packing of 
nonspherical particles, as well as flexible particles m 
and polydisperse samples m- Simulations and experi¬ 
mental work have focused on particular shapes such as 
ellipsoids mm, spherocylinders [22] [23] , and Platonic 
polyhedra [24| [2g . Convex shapes appear to achieve ran¬ 
dom packing densities well above that of spheres. This 
situation is similar to the situation for the optimal (non- 
random) packing of nonspherical particles: a conjecture, 
attributed to Ulam, posits that spheres have the lowest 
optimal packing density among all convex shapes [26] . 
Simulations with a large collection of shapes have ver¬ 
ified in each case that the shape can indeed be packed 
more densely than spheres m, and it is proven that 
the same is true for any sufficiently spherical shape, 
where asphericity is measured by 7 , the radius ratio of 


smallest circumscribed sphere to largest inscribed sphere 
[28l [29]. It has also been conjectured that all convex 
shapes below a certain asphericity (Jiao and Torquato 
suggest 7 < 1 . 2 ) have random packing densities above 
that of spheres [24l |30] . A mean-field calculation, focus¬ 
ing on axis-symmetric shapes, has provided some theo¬ 
retical backing to the conjecture isniisi]. In contrast to 
Ulam’s conjecture, the sphere cannot be a global mini¬ 
mum, since very elongated particles have random pack¬ 
ing densities in inverse proportion to their aspect ratio, 
so they can be arbitrarily small [32] . 

Moreover, unlike the optimal density, which is rig¬ 
orously defined, the random-close-packing density does 
not have an accepted mathematical definition. Oper¬ 
ationally, random close packing is often defined as the 
packing obtained through some fixed numerical proto¬ 
col [33], for example, a molecular dynamics simulation 
of slowly expanding hard particles [34] or a sequence of 
energy minimization of a system of soft particles to con¬ 
verge to the density where the repulsion energy is on the 
verge of becoming nonzero [35]. We present a heuristic 
calculation that predicts, irrespective of the preparation 
protocol, that any sufficiently spherical convex shape will 
yield a higher density than spheres would given the same 
protocol. An important assumption we make about the 
protocol is that it produces isostatic packings when ap¬ 
plied to spherical particles; that is, the number of in¬ 
terparticle contacts is equal to the number of degrees of 
freedom, as predicted by Maxwellian constraint counting 
(when unconstrained “rattlers” are removed). Isostatic- 
ity is observed in nearly all random packing protocols 
for spheres, but random packings of nonspherical parti¬ 
cles are usually hypostatic, with far fewer contacts than 
degrees of freedom [ 21 ]. 

The calculation we present is an extension to general 
convex shapes in three dimensions of the calculation per¬ 
formed by Donev et al. in two dimensions for ellipses m- 
The calculation relies on the assumption that for nearly 
spherical particles, the following two procedures would 
produce similar results: ( 1 ) applying the random pack¬ 
ing protocol to spheres and then deforming the spheres 
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to the shape of interest while maintaining a compressive 
pressure; and (2) applying the random packing proto¬ 
col directly to the particles of interest. The rationale is, 
that nearly spherical particles will behave like sphere for 
the majority of the compression protocol, and their as- 
phericity will only come into effect near the end, when 
the configuration is similar to a random packed configu¬ 
ration of spheres. This argument has been shown to hold 
true in the case of ellipses [21] . 

Therefore, we begin with a random packed configura¬ 
tion of N spheres of equal diameter a, and we continu¬ 
ously deform it into a packing of particles all congruent 
to a given body. Geometrically, the body is described as 
a compact subset K We can also define its radial 

function rx(u) = msiKxueK ^ as the distance from the 
center of K to its boundary along a certain direction u. 
Without loss of generality, we will assume that the largest 
ball contained in K is of diameter cr, and let a + Act be 
the diameter of the smallest concentric ball that contains 
K. Thus, e = Aaja = 7 — 1 measures the eccentricity of 
K. We consider a one-parameter family of nested convex 
shapes Kf, 0 < t < 1 such that Kq is the ball of diameter 
a and Ki = K. As a first step, we calculate the change 
in volume if the orientation of each particle is fixed: we 
assign to each sphere a rotation matrix Ri, i = 1,..., A, 
so that the packing at time t is given by the collection 
{RiKf + Xi(t) : i = 1,..., N}. Since the sphere is invari¬ 
ant under rotations, any choice of the matrices Ri gives 
the same configuration at t = 0. 

Throughout the deformation, we maintain a compres¬ 
sive pressure on the sample by applying external forces 
to the particles along the boundary. We label by the 
external force applied to particle i (F^ = 0 for particles 
away from the boundary). By definition, we measure the 
change in volume by the work performed by the particles 
against these forces: pAV = 

From the balance of forces on each particle at t = 0 we 
have that 

Fi = fijUij, (1) 

jedi 

where di is the set of indices of particles in contact with 
particle i at time t = 0, Uij is a unit normal vector at the 
contact point between particles i and j, and fij = fji 
is the contact force magnitude. The pairs of particles 
that are in contact, their contact normal, and contact 
force will evolve throughout the deformation, but we use 
the symbols 9 , and fij without a time argument 
to denote their values at t = 0 . We also use 
denote a sum over unordered pairs {i, j} such that j G di. 
Summing over all particles, we obtain 

pAV = '^fij{nij,Axi - Axj). ( 2 ) 

ir^j 

The pressure p is an arbitrary conversion factor between 
energy and volume dimensions, which we fix by applying 


|(2)| to an infinitesimal uniform expansion, Ax^ = ax^ and 
AV = 3 aV, to obtain 3 pV = a fij^ equivalently, 

(f) = __ ( 3 ) 

Na{\di\)i^ 

where {\di\)i is the mean coordination number, which is 
6 for isostatic packings. 

We first consider pairs of particles that are at contact 
at time t = 0 and remain so throughout the deformation. 
Their relative displacement satisfies 

Axj - Axj = rK{Ri^Uij)uij - rK{RpUji)uji - crriy, 

( 4 ) 

where Uij (w^) is a unit vector in the direction from 
the center of the particle i (j) to the point of contact 
at time t = 1 . We assume the change in the packing 
geometry is continuous, so that Uij and approach 
Auij in the limit that e ^ 0 . Therefore, as a zeroth- 
order approximation, we plug in Uij = —Uji = Uij to 
get Axi - Axj = [rxiRpnij) + rxl-Rpnij) - a] ntj. 
We estimate the error between the zeroth order approx¬ 
imation and the actual displacement by assuming self- 
consistently that the size of the displacement does not 
exceed the order 0(<je) predicted by the approxima¬ 
tion, and therefore, that Uij — Uij = 0 (e) as well. To¬ 
gether with the fact that the function rx(u) is 0(cre^/^)- 
Lipschitz continuous [ 36 ], we obtain 

(nij,Axi - Axj) = 

Ar(Rf^nij) + Ar(-Rf^nij) + 0(cre^/^), 

where Ar(u) = rx(u) — {al 2 ). 

For contacts that are broken during the deformation, 
the right hand side of |( 5 )| serves only as a lower bound, 
but we assume that the displacement is still of order 
O(cre). Since the initial number of contacts is the min¬ 
imal number required to maintain stability, the number 
of contacts broken is at most equal to the number of 
new contacts made [21]. A rough upper bound for the 
number of new contacts made is the number of pairs of 
spheres, not initially in contact, which would come to 
overlap if we dilated each sphere by a factor 1 + e. The 
number of such pairs asymptotically for small e is known 
in random packing configurations of spheres to approach 
[ 37 ] • Assuming that the forces associated with 
broken contacts are at most of average magnitude (typi¬ 
cally contacts with smaller forces will break first), we see 
that the contribution to AV^ from broken contacts in |( 5 )| 
is of order at most . Thus, to leading order in e, 

pAV = 5 ] 5 ]/,,Ar(i? 7 n,,) + 0 (ye 3 / 2 ). (6) 

i j^di 

We now consider what the change in volume would be 
if instead of holding the particle orientations fixed, we 
allowed them to rotate freely. As the system is under a 
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compressive pressure, it will tend to adopt the orienta¬ 
tions that allow it to minimize the volume. Therefore we 
assume that the change in volume will be the same as 
if we held the particle orientations fixed during the de¬ 
formation but rotated them in advance (while they are 
spheres, and their orientations is physically irrelevant) to 
the orientations that will yield the lowest volume at the 
end of the deformation. Namely, we assume that p/SV is 
equal to the minimum of the right hand side of | ( 6 ) [ over 
all choices of the rotation matrices Ri. The error term 
can be moved outside of the minimization, and we are 
left with 




min 

Rieso{^) 

for all i 


EE fijAr{Rpnij) 

i jEdi 


+ 0{Ve^^^). 


(7) 

Remarkably, each term in the outer sum being minimized 
depends only on the orientation of a single particle and 
can be minimized independently of all the other terms. 
Also, as Ri is a dummy variable, we can replace it with 
its inverse, as we do for neatness henceforth. 

The density of the final packing is given by Nvk /R(1), 
where vk is the volume of a single particle and is given 
by 


vk= rKi-afci^n 

Jues^ ( 8 ) 

= ^(cr/2)^ + 47r(cr/2)^Ar + O(cr^e^), 

where ” represents the average over the unit sphere S^. 
Rewriting m in terms of the intensive density instead 
of the extensive volume and eliminating the arbitrary 
pressure by using m we get 


A(j) 

3(po 


Ar(u) 
a 12 



(a/ 2 )(|az|),(/) 


^+0(e3/2), (g) 


where the average (•)^ represents an average over the par¬ 
ticles in the initial packing. The calculation is valid also 
in dimensions d 7 ^ 3, in which case d should be substi¬ 
tuted for 3 in |(9)[ 

If, instead of minimizing over R in |(9)[ we averaged 
over R, the first two terms on the right hand side would 
cancel. Since the mean serves an upper bound for 
the minimum, this argument immediately gives A 0 > 
-0(e3/2). is almost the claim we wish to make, 

namely, that the random packing density of any nearly 
spherical convex shape is larger than that of spheres. To 
strictly bound A(/) above 0 , we have to find a better lower 
bound on the gap between the minimum and the average. 

Consider a particle i, then we are interested in the 
minimum of gi{R) = fijAr{Rnij). We can aver¬ 

age gi{R) over all rotations that map a fixed point, say 
the north pole z, on 5'^ to a given point v to obtain a 


function over 



where R^ is the rotation mapping z to v around the axis 
perpendicular to both, and Rz{0) is the rotation about 
the 2 ;-axis by an angle of 0. Clearly, the minimum of hi{v) 
is no smaller than the minimum of gi{R). The structure 
of the linear operator : Ar(u) 1 -^ hi{w) is that of a 
convolution with a zonal (i.e. invariant over rotations that 
fix the pole) measure over S‘^ [38] . This structure allows 
us to bound the Li norm of the deviation of hi(y) from 
its mean in terms of that deviation in Ar. This procedure 
also gives the following bound on the minimum: 

min 5 fi(R) - gi{R) < minhi(v) - hi{w) 

^ " ( 11 ) 

< Ci\Ar{u) — Ar(u)|. 


We do not give a complete description of this proce¬ 
dure, as it is essentially the same as Section 4 of [28] . 
The constant q is strictly positive whenever the spher¬ 
ical harmonic expansion of the zonal measure used in 
the convolution has no terms that vanish. Therefore, 
as this will not happen generically (not to mention to 
a fraction of particles approaching one), we can safely 
assume that the average value {ci)i is strictly positive. 
Using the bound (11) in |(9)[ we immediately get that 


A(/) > c|Ar(u) — Ar(u)| for some constant c. As the 
right hand side is zero only for spheres, we have obtained 
the result we were after. 

Having obtained the theoretical result for general 
shapes and general packing protocols, we now wish to 
calculate numerical estimates for specific shapes and pro¬ 
tocols. Consider a family of convex shapes, parame¬ 
terized by some single variable, a, that includes the 
sphere of diameter a dX a = ao- Let ra{u) be the 
radial function describing the shapes, and let p(u) = 
{dra{u)/da\^+)/{a/2). We define 7 ^ as a measure of the 
slope of the random packing density subject to a given 
protocol as a function of the parameter a: 


1 d(j){a) 
^ 30 da 


( 12 ) 


We take the derivative in the direction of positive a, as 
0(a) is usually not smooth at ao- From m we have 


ri = /5(u) 


1 

mw) 



(13) 


where {')i denotes the average over particles in random 
packing configuration of spheres produced by the same 
protocol. 

We estimate g for a few examples of shape families in 
the case of the random packing protocol of Jin and Makse 
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FIG. 1. We calculate numerical slope estimates for 12 shape 
families, representatives of which are illustrated here in the 
order listed in Table I (from top, left to right). 


[39] . To numerically calculate |( 13) [ we use a dataset that 
includes the coordinates and forces of a configuration of 
10^ spheres in a periodic cubic box [39] . 

Consider first the family of prolate ellipsoid, param¬ 
eterized by aspect ratio. The infinitesimal deformation 
is p(6>, ip) = cos^ d and we can solve the minimization 
analytically: 


where Fi = ^ ^ij ^ symmetric ten¬ 

sor describing the stress on particle and Amin is its 
smallest eigenvalue. We calculate a numerical value of 
^proi = 0.202 ± 0.001. The error estimate quoted in¬ 
cludes only statistical error. There could be systematic 
error from the fact that the dataset we use is not precisely 
at the isostatic point and from finite-size effects. 

For an oblate ellipsoid, p = 1 — cos^ 6>, and so 

^ {\di\)i{f) “ 3 ' 


Numerically, we get 77obi = 0.242 ± 0.003. In the case of 
triaxial ellipsoids, we consider for each 0 < p < 1 the 
family where the principal axes are given by cr < < 

aa. Then p(6>, (j)) = cos^ 0 F p sin^ 0 cos^ 0, and we get 

'^triax ~ (f /^)'^prol “1“ A^'^obl- 

Leaving aside ellipsoids, where minimization can be 
done analytically, we consider other shape families of in¬ 
terest, illustrated in |FigureT] and defined in the appendix. 
We numerically find Ri that minimizes 
for each of the 10^ spheres in the dataset and for each 
shape family. The resulting numerical values for the slope 
p = {l/3)d(l)/da are given in Table I| Since the value of 
p depends on the way the family of shapes is parameter¬ 
ized, we also give the normalized value p/\p — p\ [40]. 


shape family 

V 

v/(\p- 

-p\) 

oblate ellipsoid 

0.242 

± 

0.003 

0.94 

± 

0.01 

prolate ellipsoid 

0.202 

± 

0.001 

0.79 

± 

0.01 

lens 

0.216 


0.004 

0.86 

± 

0.01 

spherocylinder 

0.271 

± 

0.002 

1.08 

± 

0.01 

spherodisk 

0.246 

± 

0.003 

1.36 

± 

0.02 

spindle 

0.139 

± 

0.005 

0.77 

± 

0.03 

rounded tetrahedron 

0.189 

± 

0.003 

1.45 

± 

0.02 

tetrahedral puff 

0.138 

± 

0.005 

1.06 

± 

0.04 

rounded triangle 

0.306 

± 

0.002 

1.31 

± 

0.01 

triangular spindle 

0.235 

± 

0.004 

1.01 

± 

0.02 

cubic superball 

0.142 

± 

0.001 

1.32 

± 

0.01 

octahedral superball 

0.130 

± 

0.002 

1.20 

± 

0.02 


TABLE 1. Slopes of the random close packing density of one- 
parameter families of shapes at the point corresponding to the 
sphere. The slope parameter is defined as p = {l/3)d(l)/da 
and dep ends on the parameterization. The normalized ver¬ 
sion, p/{\p — p\), where p(u) = dra{u)/da, is intrinsic to the 
shape family 



FIG. 2. Random packing density of prolate and oblate ellip¬ 
soids according to simulation results of Ref. [2T] together with 
lines illustrating predicted slopes at a = 1. Inset: Random 
packing density of octahedral and cubic superballs according 
to simulation results of Ref. m together with lines illustrat¬ 
ing predicted slopes at p = 1. 


Rounded tetrahedra appear to give the largest nor¬ 
malized improvement in packing density out of all the 
families considered. This observation is in harmony with 
the fact that regular tetrahedra, out of all convex shapes 
that have been studied, seem to have the largest random 
close packing density [42] , 

We compare our predicted slopes to simulation data for 
prolate and oblate ellipsoids [21] and for superballs [41] 
in [Figure 2] The protocols used in these simulations are 
differenent than the protocol used to obtain the data on 
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which we based our numerical calculation, so this com¬ 
parison should be interpreted with caution. It is hard 
to tell from the limited data if our calculation system¬ 
atically overestimates the slope of the density curve, or 
whether nonlinearities of the curve quickly cause the data 
to diverge from the linear estimate. To resolve this un¬ 
certainty, data for shapes in these families closer to the 
spheres will be needed. 

The method we present allows us to extend the ex¬ 
isting robust knowledge about the random packing be¬ 
havior of spheres to all sufficiently spherical shapes. We 
perform some numerical calculations for specific one- and 
two-parameter families of three-dimensional shapes, but 
our method is applicable to any nearly spherical shape in 
any number of dimensions. As such, it provides a valu¬ 
able tool to test general ideas about the random packing 
behavior of nonspherical shapes, such as the conjecture 
that spheres minimize this density among all sufficiently 
spherical convex shapes. Our calculation predicts that 
indeed spheres are a local minimum, providing substan¬ 
tial backing to the conjecture. It is worth noting that this 
claim holds irrespective of the number of dimensions. In 
this aspect, the optimal (nonrandom) packing problem is 
completely different: the d-dimensional sphere in any di¬ 
mension other than 3 does not appear to be even a local 
minimum of the optimal packing density [28] . 
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Appendix: definition of shape families 

In this appendix, we give definitions of the families of 
shapes for which the slope of the random packing density 
curve was calculated in ITable II 

An oblate ellipsoid of aspect ratio a is the region 
{{x^y^z) G : (x/a)‘^ + {v/glY + < 1}, where 

a > 1. A prolate ellipsoid is the region {{x^y^z) G : 
x‘^ +y‘^{z/a)‘^ < 1}, where a > 1. A lens is the intersec¬ 
tion of two equal-sized spheres, and a spherocylinder is 
the convex hull of two equal-sized spheres. A spherodisk 
is the convex hull of all equal-sized spheres with centers 
on a given circle, and a spindle is the intersection of such 
a family of spheres. All shapes in these six families are 
axis symmetric. 

We consider a few more bodies without axial symme¬ 
try. A tetrahedral puff is the intersection of four equal¬ 
sized spheres with centers at the corners of a regular 
tetrahedron [43]. A rounded tetrahedron is the convex 
hull of four such spheres. Similarly, a triangular spindle 
and a rounded triangle, are, respectively, the intersection 
and convex hull of three equal-sized spheres at the cor¬ 


ners of an equilateral triangle. A superball is the region of 
space determined by the inequality + < 1 

m- When p = 1, we recover the Euclidean ball. When 
1 < p < oc, we call the superballs cubic, since they 
interpolate between the ball and the cube. Similarly, su¬ 
perballs with I < p < 1 are called octahedral, as they 
interpolate between the ball and the octahedron. 

The infinitesimal deformations associated with these 
families are 

Pobl-ell(u.) = Pprol-ell (u.) = 

Plens(u) ~ Psph-cyl(u) ~ \'^z\ 

Psph-disk(u) = Pspindle(u.) = \/l 

Prnd-tet(u) = -ptet-puff(u) = . JUaX ^ U • vf ^ 

Prnd-tri(u) = “ptri-spind (u) = maX U • vf 

z—1 ^ 2 ^ o 

Pcub-sup(u) ~ Poct-sup(u) ~ ^ ^ 

i=x,y,z 

where u = {ux^Uy^ Uz) is a point on the unit sphere, 
i = 1, 2,3,4, are the four vertices of a regular tetrahedron 
inscribed in the unit sphere, and vT^, i = 1,2,3, are the 
three vertices of an equilateral triangle inscribed in the 
equator. 





